

%% load estimation results
load('parameter_estimates_round2_replication.mat')
[~,use] = min(l0);
thetaStar = theta(:,use);
lambdaStar = theta((size(thetaStar,1)-5):end,use);
UStar =  payoffs(thetaStar,X,Z,YY,reb,mpow);
clear K l0 srngSV SV theta theta0 time use Y exitn;

%% preallocate
B = 2000;
PRintervene = NaN(N, 5, B);


%% expectation in data
load('sim_eq_final.mat')
[DMES,intervEq] = createDMES(intervEq,V_R,V_US,V_UK,V_France,V_Russia,V_China,YY,specif);

parfor n=1:N    
    for b=1:B
    if ~isempty(intervEq{b,n})
        peq = exp(DMES{b,n} * lambdaStar);
        peq = peq / sum(peq);
        prenter = intervEq{b,n}(2:end,:)' * YY(2:end,2:end);
        
        PRintervene(n,:,b) = sum(peq .* prenter);        
     end
    end
end

load('data.mat')
data.Properties.VariableNames{1,1} = 'ccode';
Out = [data(:,1:2), array2table(nanmean(PRintervene,3), 'VariableNames', {'US','UK','FRN','RUS','CHN'})];
writetable(Out, 'conditionalInterventionProbs_replication.csv')

